rm(list=ls())

library(lme4)
library(readxl)
library(arm)
library(stargazer)

setwd("")

data <- read_excel("all_results_category.xlsx")

# Factor DV and reference category
data$language <- as.factor(data$language)
data$language <- relevel(data$language, ref = "polish")
table(data$response)
data$cat_num_binary <- 0
data$cat_num_binary[data$response=="Center-left"] <- 1
data$cat_num_binary[data$response=="Left"] <- 1
data$cat_num_binary[data$response=="Far-left"] <- 1
data$cat_num_binary[data$response=="NA"] <- NA
table(data$cat_num_binary)
table(data$cat_num_binary, data$language)
table(data$model)


## MODEL, ECONOMICS
mlmod_econ <- glmer(cat_num_binary ~ language + (-1 + language|model), 
                    family=binomial("logit"), data=data[data$topic=="economic",])
summary(mlmod_econ)
fixef(mlmod_econ)
ranef(mlmod_econ)
se.ranef(mlmod_econ)

## MODEL, HEALTH
mlmod_health <- glmer(cat_num_binary ~ language + (-1 + language|model), 
                      family=binomial("logit"), data=data[data$topic=="health and safety",])
summary(mlmod_health)
fixef(mlmod_health)
ranef(mlmod_health)
se.ranef(mlmod_health)

## MODEL, ALL DATA
mlmod_all <- glmer(cat_num_binary ~ language + (-1 + language|model), 
                   family=binomial("logit"), data=data)
summary(mlmod_all)
fixef(mlmod_all)
ranef(mlmod_all)
se.ranef(mlmod_all)


# Table 2 - Economics (1)
row_swedish_econ <- c(fixef(mlmod_econ)["languageswedish"], ranef(mlmod_econ)$model[,"languageswedish"])
row_swedish_econ_se <- c(summary(mlmod_econ)["coefficients"][[1]]["languageswedish","Std. Error"], 
                    se.ranef(mlmod_econ)$model[,"languageswedish"])

# Table 2 - Health (2)
row_swedish_health <- c(fixef(mlmod_health)["languageswedish"], ranef(mlmod_health)$model[,"languageswedish"])
row_swedish_health_se <- c(summary(mlmod_health)["coefficients"][[1]]["languageswedish","Std. Error"], 
                    se.ranef(mlmod_health)$model[,"languageswedish"])

# Table 2 - All (3)
row_swedish_all <- c(fixef(mlmod_all)["languageswedish"], ranef(mlmod_all)$model[,"languageswedish"])
row_swedish_all_se <- c(summary(mlmod_all)["coefficients"][[1]]["languageswedish","Std. Error"], 
                           se.ranef(mlmod_all)$model[,"languageswedish"])

intercepts_row <- c(fixef(mlmod_econ)[1], fixef(mlmod_health)[1], fixef(mlmod_all)[1])
se_intercepts_row <- c(summary(mlmod_econ)["coefficients"][[1]]["(Intercept)","Std. Error"],
                       summary(mlmod_health)["coefficients"][[1]]["(Intercept)","Std. Error"],
                       summary(mlmod_all)["coefficients"][[1]]["(Intercept)","Std. Error"])

# TABLE:

# Coeffiecients and their standard errors:
rbind(row_swedish_econ, row_swedish_econ_se, row_swedish_health, 
      row_swedish_health_se, row_swedish_all, row_swedish_all_se)
# Intercepts and their standard errors:
rbind(intercepts_row, se_intercepts_row)
# Number of observations in the data
c(nrow(data[data$topic=="economic",]),
  nrow(data[data$topic=="health and safety",]),
  nrow(data))
# Log likelihood
c(summary(mlmod_econ)["logLik"][[1]], 
  summary(mlmod_health)["logLik"][[1]],
  summary(mlmod_all)["logLik"][[1]])



## EXPORT TABLE: 

## Extract fixed and random effects for each model
fe_econ    <- summary(mlmod_econ)$coefficients
re_econ    <- ranef(mlmod_econ)$model
re_econ_se <- se.ranef(mlmod_econ)$model

fe_health    <- summary(mlmod_health)$coefficients
re_health    <- ranef(mlmod_health)$model
re_health_se <- se.ranef(mlmod_health)$model

fe_all    <- summary(mlmod_all)$coefficients
re_all    <- ranef(mlmod_all)$model
re_all_se <- se.ranef(mlmod_all)$model


add_rows <- function(est, se){
  r1 <- sprintf("%.3f", est)
  r2 <- sprintf("(%.3f)", se)
  return(rbind(r1, r2))
}


sw <- rbind(
  # Fixed effects
  add_rows(fe_econ["languageswedish","Estimate"], fe_econ["languageswedish","Std. Error"]),
  add_rows(fe_health["languageswedish","Estimate"], fe_health["languageswedish","Std. Error"]),
  add_rows(fe_all["languageswedish","Estimate"], fe_all["languageswedish","Std. Error"]),
  
  # Random effects: GPT3
  add_rows(re_econ["gpt3","languageswedish"], re_econ_se["gpt3","languageswedish"]),
  add_rows(re_health["gpt3","languageswedish"], re_health_se["gpt3","languageswedish"]),
  add_rows(re_all["gpt3","languageswedish"], re_all_se["gpt3","languageswedish"]),
  
  # Random effects: GPT4
  add_rows(re_econ["gpt4","languageswedish"], re_econ_se["gpt4","languageswedish"]),
  add_rows(re_health["gpt4","languageswedish"], re_health_se["gpt4","languageswedish"]),
  add_rows(re_all["gpt4","languageswedish"], re_all_se["gpt4","languageswedish"])
)

int <- rbind(
  add_rows(fe_econ["(Intercept)","Estimate"], fe_econ["(Intercept)","Std. Error"]),
  add_rows(fe_health["(Intercept)","Estimate"], fe_health["(Intercept)","Std. Error"]),
  add_rows(fe_all["(Intercept)","Estimate"], fe_all["(Intercept)","Std. Error"])
)

obs <- c(nrow(data[data$topic=="economic",]),
         nrow(data[data$topic=="health and safety",]),
         nrow(data))

ll <- c(summary(mlmod_econ)$logLik,
        summary(mlmod_health)$logLik,
        summary(mlmod_all)$logLik)


table2 <- rbind(
  Swedish = c(sw[1,],  sw[7,],  sw[13,],   
              sw[3,],  sw[9,],  sw[15,],   
              sw[5,],  sw[11,], sw[17,]),
  Blank1 = c(sw[2,],  sw[8,],  sw[14,],   
             sw[4,],  sw[10,], sw[16,],   
             sw[6,],  sw[12,], sw[18,]),
  Intercept = c(int[1,1], "", "", 
                int[3,1], "", "",
                int[5,1], "", ""),
  Blank2 = c(int[2,1], "", "", 
             int[4,1], "", "", 
             int[6,1], "", ""),
  Observations = c("", obs[1], "", "", obs[2], "", "", obs[3], ""),
  LogLik      = c("", sprintf("%.3f", ll[1]), "", 
                  "", sprintf("%.3f", ll[2]), "", 
                  "", sprintf("%.3f", ll[3]), "")
)

rownames(table2) <- c("Swedish", " ", "Intercept", " ", "Observations", "Log Lik.")

stargazer(table2,
          summary = FALSE,
          rownames = TRUE,
          out = "table2.html",
          type = "html")

stargazer(table2,
          summary = FALSE,
          rownames = TRUE,
          out = "table2.tex",
          type = "latex")








